import numpy
import Image as pil
#import pylab



def load(filename):
    """ pil load """
    i=pil.open(filename)
    ri=numpy.asarray(i)
    return ri

def resize(i, s):
    """ pil resize """
    im1= pil.fromarray(numpy.uint8(i))

    im2= im1.resize(s, pil.NEAREST)
    ri= numpy.asarray(im2)
    return ri

def convert_to_grayscale(im):
    """ convert a rgb-image to a grayscale image. it is done by addition of the three channels."""
    im = im.astype(float)
    s = im.shape
    im2 = numpy.zeros([s[0], s[1]])
    #for x in range(0, s[0]):
    #    for y in range(0, s[1]):
    #        im2[x, y] = sum(im[x, y])
    if len(s) >= 3:
        im2[:,:] = im[:,:, 0] + im[:,:, 1] + im[:,:, 2]
    else:
        im2=im
    im2 = im2 / (3*255.0)
    return im2

def is_grayscale(im):
    """ checks wether image is a grayscale or rgb"""
    if len(im.shape)>=3:
        return 0
    else:
        return 1

def find_horizontal_edges(im):
    """simply comutes the edges by difference of the rows"""
    im = im.astype(float)

    s = im.shape
    im2 = numpy.zeros(s)
    if is_grayscale(im):
        im3 = numpy.zeros(s)
        im3[:, 1:]=im[:,0:(s[1]-1)]
        im2=im-im3
    #im2 = im2 + abs(im2.min())
    #im2 = im2 / float(im2.max())
    return im2
	
def find_vertical_edges(im):
    """simply comutes the edges by difference of the columns"""
    im = im.astype(float)

    s = im.shape
    im2 = numpy.zeros(s)
    if is_grayscale(im):
        im3 = numpy.zeros(s)
        im3[1:,:] = im[0:(s[0]-1), :]
        im2 = im-im3
    #im2 = im2 + abs(im2.min())
    #im2 = im2 / float(im2.max())
    return im2

def histogramm(im):
    maxi=im.max()
    if maxi > 0:
        i2=im/float(maxi)*255
    else:
        i2=im
    i2=i2.astype(int)
    h = numpy.zeros(256)
    if is_grayscale(im):
        s = im.shape
        for x in xrange(0, s[0]):
            for y in xrange(0, s[1]):
                h[int(im[x, y])] += 1
    return h

def binarize(im, threshold):
    im = im.astype(float)
    im2 = (im < threshold)
    im2 = im2.astype(float)
    return im2

def calc_row_average_horizontal(im):
	im=im.astype(float)
	s=im.shape
	im2=numpy.zeros(s[0])
	if len(s) <=2:
		for x in range(0, s[0]):
			im2[x]=sum(im[x,:])
		im2/=float(s[1])
	else:
		im2= None
	return im2

def calc_row_average_vertical(im):
	im=im.astype(float)
	s=im.shape
	im2=numpy.zeros(s[1])
	if len(s) <=2:
		for x in range(0, s[1]):
			im2[x]=sum(im[:,x])
		im2/=float(s[0])
	else:
		im2=None
	return im2

def calc_row_sum_horizontal(im):
    im=im.astype(float)
    s=im.shape
    im2=numpy.zeros(s[0])
    if len(s) <=2:
        for y in range(0, s[0]):
            im2[y]=sum(im[y,:]) 
    else:
        im2=None
    return im2

def get_histo_entropy(v):
    h=scipy.histogram(v, 256, normed=True)
    e=0
    #print h
    for n in xrange(0, len(h)):
        if h[0][n]!= 0:
            e+=(-h[0][n]*scipy.log(h[0][n]))
    return e

def get_entropy(v):
    e=0
    max_v=float(sum(v))
    for n in xrange(0, len(v)):
        p=v[n]/max_v
        if p > 0:
            e+=(-p*scipy.log(p))
    return e

def calc_correlation(x, y):
	mx=x.mean()
	my=y.mean()
	s=x.shape[0]
	c=1/float(s)*sum((x-mx)*(y-my))
	c=c/((1/float(s)*sum((x-mx)**2))**0.5*(1/float(s)*sum((y-my)**2))**0.5)
	return c
	
def thin_pic(im):
	s=im.shape
	x=0

	y=0
	
	im2=numpy.zeros(s)
	for  x in range(1, s[1]):
		for y in range(1, s[0]):
			if im[y, x]>=1:
				im3=im[y-1:y+2, x-1:x+2]
				ss=sum(sum(im3))
				#print ss
				if (ss >= 3):
					im2[y, x]=1.0
	return im2

def distance_transform(im):
    s=im.shape
    index_range=list([-1, 0, 1])
    increment=0
    ok=1
    while ok:
        ok = 0
        #print ".",
        im2 = im.copy()
        increment+=1
        for  x in xrange(1, s[1]-1):
            for y in xrange(1, s[0]-1):
                do_inc=1
                for xi in index_range:
                    for yi in index_range:
                        if im[y-yi, x-xi] != increment:
                            do_inc=0
                if do_inc:
                    im2[y, x]=increment+1


        d = abs(im-im2)
        ok = sum(sum(d))
        im = im2.copy()
        #number_of_iterations += 1
    return im2

def skeleton_pic(im):
    s = im.shape
    im4=numpy.zeros(s)
    im2=distance_transform(im)
    index_range = list([-1, 0, 1])
    for  x in xrange(1, s[1]-1):
        for y in xrange(1, s[0]-1):
            im3=im2[y-1:y+2, x-1:x+2]
            #print im3
            m=im3.max()
            if (im2[y, x]==m) and (im2[y, x] > 0):
                im4[y, x]=1
            else:
                im4[y, x]=0
    return im4

def game_of_life(im):
	s=im.shape
	im_n=numpy.zeros(s)
	for  x in range(1, s[1]):
		for y in range(1, s[0]):
			im3=im[y-1:y+2, x-1:x+2]
			#print im3.shape
			ss=sum(sum(im3))
			if (ss < 2) or (ss> 3):
				im_n[y, x]=0
			if (2 <= ss <= 3) and (im[y, x]>=1.0):
				im_n[y, x]=1
			if (ss == 3) and (im[y, x]==0):
				im_n[y, x]=1
	return im_n
	
def flood_fill(im, x, y, target, replacement):
    "Uses the flood_fill to fill image with a color defined by replacement, when pixel is of target-color."
    q = list()
    if im[y, x] != target:
        return im
    im2 = im.copy()
    if is_grayscale(im):
        #print replacement
        q.append([x, y])
        #im[y, x]=replacement
        s = im2.shape
        while len(q) > 0:
            #print q
            x, y = q[0]
            if im2[y, x] == target:
                im2[y, x] = replacement
            q = q[1:]
            if (x-1  >= 0):
                if im2[y, x-1] == target:
                    q.append([x-1, y])
                    im2[y, x-1] = replacement
            if (x + 1  < s[1]):
                if im2[y, x + 1] == target:
                    q.append([x + 1, y])
                    im2[y, x + 1] = replacement
            if (y + 1 < s[0]):
                if im2[y + 1, x] == target:
                    q.append([x, y + 1])
                    im2[y + 1, x] = replacement
            if (y-1 >= 0):
                if im2[y-1, x] == target:
                    q.append([x, y-1])
                    im2[y-1, x] = replacement
    return im2
	
def segment_picture(im, color):
    """ Selects all pixels of one object in image and uses one value per object to fill it."""
    im2 = im.copy()
    if is_grayscale(im):
        s = im.shape
        rep = color + 1
        for x in xrange(0, s[1]):
            for y in xrange(0, s[0]):
                #print y, x
                if im2[y, x] == color:
                    #print y, x
                    im2 = flood_fill(im2, x, y, color, rep)
                    #pylab.matshow(im2)
                    rep += 1
    return im2

def reduce_size_roi(im):
    """ reduce size of image by cutting background"""
    minx, miny = im.shape[1], im.shape[0]
    if (minx!=0) and (miny!=0):
        maxx, maxy = 0,0 
        minx, miny, maxx, maxy=get_koordinates_roi(im)
        #print minx, miny, maxx, maxy
        i2=im[miny:maxy+1, minx:maxx+1].copy()
    else:
        i2=im
    return i2

def get_koordinates_roi(im):
    """ reduce size of image by cutting background"""
    minx, miny = 0,0
    maxx, maxy = 0,0 
    #res=[0,0,0,0]
    a=numpy.array(im.nonzero())
    if (a.shape[0]>0) and (a.shape[1]>0):
        minx=a[1,:].min()
        miny=a[0,:].min()
        maxx=a[1,:].max()
        maxy=a[0,:].max()   
    return [minx, miny, maxx, maxy]

def white_balance(im):
    """ pixel values are divided by the maximum value to spread the values in the interval [0..1]"""

    im2=im/float(im.max())
    return im2

def maximise_contrast(im):
    """ contrast is spread in the interval [0..1]"""
    m=im.min()
    im2=im-m
    im2=im2/float(im2.max())
    return im2

def filter(mat, f):
    xsize, ysize= mat.shape
    xrange=range(1, xsize-1)
    yrange=range(1, ysize-1)
    mat2=scipy.zeros([xsize, ysize])
    mat3=scipy.zeros([3, 3])
    for x in xrange:
        for y in yrange:
            mat3=mat[x-1:x+2, y-1:y+2]
            mat3=mat3*f
            mat2[x ,y]=sum(sum(mat3))
    return mat2

def hough_transformation(im, c=0, max_alpha=100):
    """ The hough transfomration is calculated with the angel running from 0..pi. The corresponding line is of the form: d=x*cos(alpha)+y*sin(alpha)."""
    xsize, ysize=im.shape
    max_d=((xsize)**2+(ysize)**2)**0.5
    pi_vec=numpy.linspace(0, numpy.pi, max_alpha)
    akk=numpy.zeros([len(pi_vec), 2*max_d])
    for x in xrange(0, xsize):
        for y in xrange(0, ysize):
            if im[x, y] >= c:
                for ai in xrange(0, len(pi_vec)):
                    alpha=pi_vec[ai]
                    d=int(x*numpy.cos(alpha)+y*numpy.sin(alpha)+max_d)
                    #print d
                    akk[ai, d]+=1
    return akk

def double_image_size(im):
    """The image size is doubled. All missing pixels are calculated through linear interpolation."""
    im=im.astype(float)
    s=im.shape
    #print s
    im2=numpy.zeros([2*s[0], 2*s[1]])
    for x in xrange(0, s[0]-1):
        for y in xrange(0, s[1]-1):
            v=im[x, y]
            im2[2*x, 2*y]=v

            im2[2*x+1, 2*y]=(v+im[x+1,y])/2.0
            im2[2*x+1, 2*y+1]=(v+im[x+1, y+1])/2.0
            im2[2*x, 2*y+1]=(v+im[x, y+1])/2.0
    for x in xrange(0, s[0]):
        im2[2*x, 2*s[1]-2]=im[x, s[1]-1]
        im2[2*x+1, 2*s[1]-2]=im[x, s[1]-1]
        im2[2*x, 2*s[1]-1]=im[x, s[1]-1]
        im2[2*x+1, 2*s[1]-1]=im[x, s[1]-1]
    for y in xrange(0, s[1]):
        im2[2*s[0]-2, 2*y]=im[s[0]-1, y]
        im2[2*s[0]-2, 2*y+1]=im[s[0]-1, y]
        im2[2*s[0]-1, 2*y]=im[s[0]-1, y]
        im2[2*s[0]-1, 2*y+1]=im[s[0]-1, y]
    #im2[2*s[0]-1, :]=im2[2*(s[0]-1), :]
    #im2[:, 2*s[1]-1]=im2[:, 2*(s[1]-1)]
    return im2

def calc_alpha(i1, i2):
    m1=i1.mean()
    m2=i2.mean()
    #print m1, m2
    alpha=scipy.log(m1)/scipy.log(m2)
    return alpha

def mosaic(i, xc, yc):
    s=i.shape
    i2=numpy.zeros(s)
    xsize=s[1]/xc
    ysize=s[0]/yc
    #print xsize
    i3=scipy.misc.imresize(i, [ysize, xsize])
    i3=white_balance(i3)
    #print i3
    for u in xrange(0, xc):
        for v in xrange(0, yc):
            alpha=calc_alpha(i[v*ysize:(v+1)*ysize, u*xsize:(u+1)*xsize], i3)
            #print alpha
            for x in xrange(0, xsize):
                for y in xrange(0, ysize):
                  i2[v*ysize+y, u*xsize+x]=(i[v*ysize+y, u*xsize+x]+i3[y, x]**alpha)/2.0
    
    return i2

def jaehne_test_pattern(s):
    x=scipy.arange(-s, s)
    y=scipy.arange(-s, s)
    mx,my=scipy.meshgrid(x,y)
    r=(mx**2+my**2)**0.5
    km=0.8*scipy.pi
    rm=s
    w=rm/10
    t1=scipy.sin( (km * r**2) / (2 * rm) )
    t2=0.5*scipy.tanh((rm - r)/w) + 0.5
    g=(t1+t2+1)/2
    return g

def transform_full_circle(im):
    mx=im.shape[0]
    my=im.shape[1]
    ri=scipy.zeros([2*my, 2*my])
    for x in xrange(0, 2*my):
        for y in xrange(0, 2*my):
            r=((y-my)**2)**0.5
            r2=((x-my)**2+(y-my)**2)**0.5
           
            if r2 > 0:
                if (y-my) > 0:
                    phi2=scipy.arccos((x-my)/r2)
                else:
                    phi2=2*scipy.pi-scipy.arccos((x-my)/r2)
                ox=int(phi2/(2*scipy.pi)*mx)
                #print (y-my)/r2
                oy=int(r2)
            else:
                ox=0
                oy=0
            phi=x/float(mx)*scipy.pi
            x2=int(r*scipy.cos(phi)+my)
            y2=int(r*scipy.sin(phi))
            #print ox, oy, mx, my
            if ox < im.shape[0] and oy < im.shape[1]:
                #print r, phi, x2, y2,x,y
                ri[x, y]=im[ox,oy]
    return ri

def transform_half_circle(im):
    mx=im.shape[0]
    my=im.shape[1]
    ri=scipy.zeros([2*my, my])
    for x in xrange(0, 2*my):
        for y in xrange(0, my):
            r=((y-my)**2)**0.5
            r2=((x-my)**2+(y)**2)**0.5
           
            if r2 > 0:
                

                phi2=scipy.arccos((my-x)/r2)

                ox=int(phi2/(scipy.pi)*mx)
                #print (y-my)/r2
                oy=int(r2)
            else:
                ox=0
                oy=0
            phi=x/float(mx)*scipy.pi
            x2=int(r*scipy.cos(phi)+my)
            y2=int(r*scipy.sin(phi))
            #print ox, oy, mx, my
            if ox < im.shape[0] and oy < im.shape[1]:
                #print r, phi, x2, y2,x,y
                ri[x, y]=im[ox,oy]
    return ri

def smooth_image(i):
    f=1/9.0*scipy.ones([3,3])
    r=scipy.ndimage.convolve(i, f)
    return r

def lighten_image(i):
    f=scipy.zeros([3,3])
    f[1,1]=1.2
    r=scipy.ndimage.convolve(i, f)
    return r

def sharpen_image(i): 
    f=-1*scipy.ones([3,3])
    f[0,0]=0
    f[2,0]=0
    f[0,2]=0
    f[2,2]=0  
    f[1,1]=5
    #print f
    r=scipy.ndimage.convolve(i, f)
    return r

def fisheye(im,mx, my,rg):
    #my=im.shape[0]/2
    #mx=im.shape[1]/2
    ri=scipy.zeros([im.shape[0], im.shape[1]])
    for x in xrange(0, im.shape[1]):
        for y in xrange(0, im.shape[0]):
            r=((y-my)**2+(x-mx)**2)**0.5
            
           
            if (x-mx)**2+(y-my)**2 < rg**2:
                alpha=scipy.arcsin((y-my)/(r+1))
                beta=scipy.arccos((x-mx)/(r+1))
                gamma=scipy.arcsin(r/rg)
                r2=2*gamma/scipy.pi*rg
                x2=int(r2*scipy.cos(beta)+mx)
                y2=int(r2*scipy.sin(alpha)+my)
                #print ox, oy, mx, my
                if 0<x2 < im.shape[1] and 0<y2 < im.shape[0]:
                    #print r, phi, x2, y2,x,y
                    ri[y, x]=im[y2,x2]

            else:
                ri[y,x]=im[y,x]
    return ri

def miniturize(i):
    mat=scipy.ones([3,3])
    
    ydim=i.shape[0]
    xdim=i.shape[1]
    mx=xdim/2
    my=ydim/2
    ri=scipy.zeros([ydim, xdim])
    
    for y in xrange(0, ydim):
        r=abs(y-my)
        filter_size=int(8*r/my)+1
        mat=scipy.ones([filter_size,filter_size])/float(filter_size**2)
        #print mat.shape
        filtered_image=scipy.ndimage.convolve(i, mat)
        #print filtered_image.shape, i.shape
        for x in xrange(0, xdim):
            ri[y,x]=filtered_image[y,x]

    return ri

    